% This script is responsible for the calculating E-k relations for all
% valence subbands. The main outline of the script is:
% 1) LuttingerParameters.m is run to obtain the Luttinger parameters.
% 2) The zone center (k=0) subband energies are determined by searching
%    an energy grid ranging over the energies allowing a bound solution.
%    The solutions are refined up to a certain energy accuracy E_acc.
% 3) Two plots of the valence band structure are shown, distiguishing
%    between heavy hole and light hole components of the wavefunctions.
% 4) The wavefunction data is saved as an ascii file under the name
%    valence_wf_<subband number>.dat.
% 5) The dispersion relation is found by starting from the zone center
%    energy, and increasing k to find a new value of E(k). The previous
%    values for E(k) are then used to predict a value which is used as an
%    initial guess in the next k step.
% 6) The wavefunctions for the solutions are calculated and used to find
%    the relative interband optical transition strenght between each
%    calculated valence subband and a specified number of conduction subbands.
%    The relative strength data is saved as mat files.

clear all; close all; clc;

% Define global variable
global z  hbar  meV  V  dz  N  gamma1  gamma2  gamma3;

%% General definitions

hbar      = 1.05e-34;          % J*sec
meV       = 1.602e-22;         % K
m0        = 9.1e-31;           % free e mass (kg)
hbar      = 1.05e-34;          % h-bar (J*sec)

%% Reading the Al content in AlGaAs - generated by PH codes

xr    = load('x.r');        
z     = xr(:,1);               % z grid (the growth direction) (m^-10)
x     = xr(:,2)';               % Al alloy concentration in GaAs/AlGaAs HS

%% Masses

m_hh     = 0.5 + 0.29*x;       % heavy holes mass (kg)
m_lh     = 0.087 + 0.063*x;    % heavy holes mass (kg)
m_hh_110 = 0.85*m_hh./m_hh;    % hh mass along <110> direction (kg)

%% Luttinger parameters for AlGaAs

coeff    = [1 -2  0;
            1  2  0;
            1  0 -2];
const    = hbar^2/(2*m0) * [1./m_hh; 1./m_lh; 1./m_hh_110];                 
gamma    = (coeff^(-1) * const)';

gamma1   = gamma(:,1);
gamma2   = gamma(:,2);
gamma3   = gamma(:,3);

% Potential profile (from the PH codes)
Vr        = load('v.r');
V         = Vr(:,2);

% % Luttinger params
% LuttingerParams();
% 
% gamma1r   = load('gamma1.r');
% gamma1    = gamma1r(:,2);
% gamma2r   = load('gamma2.r');
% gamma2    = gamma2r(:,2);
% gamma3r   = load('gamma3.r');
% gamma3    = gamma3r(:,2);

%% Simulation parameters

dz                = z(2)-z(1);      % grid step size
N                 = size(V,1);   	% num of grid points
nrm               = 1e7;

min_E_diff        = 0.05*meV;       % helps to determine whether two solutions are the same
a_lim             = 1e5;            % max. for a
k_t               = 1e5;            % k_t for band-edge calculation
E_min             = min(V);         % energy minimum
E_max             = min(V(1),V(N)); % energy maximum
E_step            = 0.1*meV;    	% energy step for band-edge calculation
E_acc             = meV*1e-9;       % energy resolution of the calculated subband energies

% Search parameters
subband_E_search  = 0.15*meV;       % energy search range
E_abort           = meV;            % abort search when E range > E_abort
k_t_max           = 1e9;            % max. calculated k_t
k_t_step          = 2e6;            % k_t step for E-k dispersion
num_cb_subbands   = 4;              % number of conduction band subbands taken into account
%num_vb_subbands   = nr-1;        % number of valence band subbands taken into account

% Plotting parameters
cl = ['m  ';'c  ';'r  ';'g  ';'b  ';'y  ';'m-.';'c-.';'r-.';'g-.';'b-.';'y-.';'m  ';'c  ';'r  ';'g  ';'b  ';'y  ';'m-.';'c-.';'r-.';'g-.';'b-.';'y-.';'m  ';'c  ';'r  ';'g  ';'b  ';'y  '];
ch = cl;

%% Calculating the band-edge subband energies and wavefunctions

disp('Calculating the band-edge subband energies and wavefunctions:');

% Plot the potential
fig_l = figure;
plot(z/1e-10, V/meV, 'k');
hold; title('LH Subbands');
fig_h = figure;
plot(z/1e-10, V/meV, 'k');
hold; title('HH Subbands');

nr           = 1;
old_sign_l   = 1;
old_sign_h   = 1;
old_sign_d   = 1;
old_sign_d2  = 1;
old_sign_d3  = 1;
ii           = 1;

% Staring energy loop
for (E = E_min:E_step:E_max)
    disp(['ii=' num2str(ii) ', E=' num2str(E)]);

    % The transfer matrix
    tf_mat           = TransferMatrix(E, k_t);
    h_10(ii)         = tf_mat(3,3);
    l_10(ii)         = tf_mat(4,3);
    h_01(ii)         = tf_mat(3,4);
    l_01(ii)         = tf_mat(4,4);
    c_min(ii)        = -(tf_mat(1,3)*tf_mat(1,4) + tf_mat(2,3)*tf_mat(2,4)) / (tf_mat(1,4)^2 + tf_mat(2,4)^2);
    diff_h_min(ii)   = tf_mat(3,3) + c_min(ii)*tf_mat(3,4);
    diff_l_min(ii)   = tf_mat(4,3) + c_min(ii)*tf_mat(4,4);
    d_min(ii)        = diff_h_min(ii)^2 + diff_l_min(ii)^2;
    dum(ii)          = abs(h_10(ii)*h_01(ii) + l_01(ii)*l_10(ii));

    if (ii>1)
        diff_d_min(ii)  = d_min(ii) - d_min(ii-1);
        new_sign_d      = sign(diff_d_min(ii));

        if (ii>2)
            if (new_sign_d-old_sign_d==2)

                [a, E_result] = WavefunctionMinimum(E-2*E_step,E,k_t,E_acc);

                disp(['E_result=' num2str(E_result/meV)]);

                if (nr==1) | ( (nr>1) & abs(E_result-E_state(nr-1))>min_E_diff )
                    E_state(nr) = E_result;
                    [f_h, f_l]  = PlotWavefunction(E_result,k_t,1,a);

                    %                     norm = sqrt(trapz(z, f_l.^2+f_h.^2));
                    %                     f_h = f_h ./ norm;
                    %                     f_l = f_l ./ norm;

                    disp(['Edge=' num2str((f_l(N)^2 + f_h(N)^2)/(max(abs(f_l))^2+max(abs(f_h))^2))]);

                    if ( (f_l(N)^2 + f_h(N)^2)/(max(abs(f_l))^2+max(abs(f_h))^2) < 1e-4 )
                        E_state(nr) = E_result;
                        coeff(nr)   = a;
                        wf_l(nr,:)  = f_l;
                        wf_h(nr,:)  = f_h;

                        % Plotting
                        figure(fig_l);
                        plot(z/1e-10, E_state(nr)/meV + (wf_l(nr,:).^2)/nrm, cl(nr,:));
                        drawnow;
                        refresh(fig_l);
                        figure(fig_h);
                        plot(z/1e-10, E_state(nr)/meV + (wf_h(nr,:).^2)/nrm, ch(nr,:));
                        drawnow;
                        refresh(fig_h);

                        % Saving to file
                        wf_save = [z, f_h', f_l'];
                        eval(['save wf_h' num2str(nr) '.dat wf_save -ascii']);

                        nr = nr + 1;
                    end
                end
            end
        end

        old_sign_d = new_sign_d;
    end

    if (ii>1)
        diff3(ii)   = dum(ii) - dum(ii-1);
        new_sign_d3 = sign(diff3(ii));
        if (ii>2)
            if (new_sign_d3-old_sign_d3==2)

                if (nr==1) | ( (nr>1) & abs(E_result-E_state(nr-1))>min_E_diff )

                    E_state(nr) = E_result;
                    [f_h,f_l]   =  PlotWavefunction(E_result,k_t,1,a);

                    disp(['Edge=' num2str((f_l(N)^2 + f_h(N)^2)/(max(abs(f_l))^2+max(abs(f_h))^2))]);

                    if ( (f_l(N)^2 + f_h(N)^2)/(max(abs(f_l))^2+max(abs(f_h))^2) < 1e-4  )
                        E_state(nr) = E_result;
                        coeff(nr)   = a;
                        wf_l(nr,:)  = Normalize(f_l);
                        wf_h(nr,:)  = Normalize(f_h);

                        % Plotting
                        figure(fig_l);
                        plot(z/1e-10, E_state(nr)/meV + (wf_l(nr,:).^2)/nrm, cl(nr,:));
                        drawnow;
                        refresh(fig_l);
                        figure(fig_h);
                        plot(z/1e-10, E_state(nr)/meV + (wf_h(nr,:).^2)/nrm, ch(nr,:));
                        drawnow;
                        refresh(fig_h);

                        % Saving to file
                        wf_save = [z, f_h', f_l'];
                        eval(['save wf_h' num2str(nr) '.dat wf_save -ascii']);

                        nr = nr + 1;

                    else

                        [a, E_result] = WavefunctionMinimum(E_result, E+2*E_step, k_t, E_acc, E_step/50, 2);
                        E_state(nr) = E_result;
                        [f_h,f_l]   =  PlotWavefunction(E_result,k_t,1,a);

                        disp(['Edge=' num2str((f_l(N)^2 + f_h(N)^2)/(max(abs(f_l))^2+max(abs(f_h))^2))]);

                        if ( (f_l(N)^2 + f_h(N)^2)/(max(abs(f_l))^2+max(abs(f_h))^2) < 1e-4 )
                            E_state(nr) = E_result;
                            wf_l(nr,:)  = Normalize(f_l);
                            wf_h(nr,:)  = Normalize(f_h);

                            % Plotting
                            figure(fig_l);
                            plot(z/1e-10, E_state(nr)/meV + (wf_l(nr,:).^2)/nrm, cl(nr,:));
                            drawnow;
                            refresh(fig_l);
                            figure(fig_h);
                            plot(z/1e-10, E_state(nr)/meV + (wf_h(nr,:).^2)/nrm, ch(nr,:));
                            drawnow;
                            refresh(fig_h);

                            % Saving to file
                            wf_save = [z, f_h', f_l'];
                            eval(['save wf_h' num2str(nr) '.dat wf_save -ascii']);

                            nr = nr + 1;
                        end
                    end
                end
            end
        end

        old_sign_d3 = new_sign_d3;
    end

    E_index(ii)  = E;
    ii           = ii+1;
end   % -- end of for loop --

%% Calculating the E-k dispersion relations

clear wav_con  rts_TE  rts_TM;

if ~exist('num_vb_subbands')
    num_vb_subbands = nr - 1;
end

% Conduction band wavefunctions
for (cb_num = 1:num_cb_subbands)
    eval(['wf=load(''wf_e' num2str(cb_num) '.r'');']);
    wavcon(:,cb_num) = wf(:,2);
end  % -- end of for loop --

disp('Calculating the E-k dispersion relations:');

% Valence band calculations
for (nr = 1:num_vb_subbands)
    clear  subband_E  coeff;

    % Definitions
    E_0               = E_state(nr);                   % subband edge energy
    subband_E(1)      = E_0;
    sbnd_E_srch_act   = subband_E_search;              % search range
    sbnd_E_srch_old   = sbnd_E_srch_act;               % old search range (previous k)
    E_step            = subband_E_search/10;           % energy step in range
    E_step_old        = E_step;                        % old energy step
    dE                = 0;                             % difference between two previous energies
    jj                = 1;
    abort             = 0;
    error             = 0;                             % 1 if no minimum found
    sbnd_E_lim        = meV*3e-6;                      % minimum search range
    old_error         = 0;                             % diff between the predicted and calculated E
    accuracy_modifier = 1;
    min_num           = 1;                             % number of minima in search range
    E_acc_act         = E_acc;
    k_t_mat           = [(0:0.25:1.75)  2:(k_t_max/k_t_step)]*k_t_step;

    tic;
    
    % Iterating over k_t
    while (jj < size(k_t_mat,2)) & ~abort

        k_t                 = k_t_mat(jj+1);

        min_E               = subband_E(jj) + dE - sbnd_E_srch_act;
        max_E               = subband_E(jj) + dE + sbnd_E_srch_act;
        min_E_arc(jj+1)     = min_E;
        max_E_arc(jj+1)     = max_E;

        [a,E_result,error]  = WavefunctionMinimum(min_E,max_E,k_t,E_acc_act,E_step, min_num);

        if (jj/50==floor(jj/50))
            drawnow;
        end

        if (error)                       % -- minimum not found

            if (min_num==2) & (jj<4)
                min_num = 1;
            else
                num_steps = 2*sbnd_E_srch_act/E_step;

                if (E_step>E_step_old/10) & (sbnd_E_srch_act<sbnd_E_srch_old*10)

                    if (num_steps<3e2)
                        E_step    = E_step/2;
                        num_steps = sbnd_E_srch_act/E_step;
                        disp(['Adjusting E_step to ', num2str(E_step/meV), ', ', num2str(floor(num_steps)), ' steps']);
                    else
                        sbnd_E_srch_act = sbnd_E_srch_act*2;
                        E_step          = sbnd_E_srch_act/10;
                        disp(['Too many steps - increaing search area to ', num2str(sbnd_E_srch_act/meV), ', E_step=', num2str(E_step/meV), ', 10 steps']);
                    end

                elseif (sbnd_E_srch_act<5*sbnd_E_srch_old)

                    sbnd_E_srch_act = sbnd_E_srch_old*5;
                    E_step          = E_step_old/2;
                    num_steps       = sbnd_E_srch_act/E_step;
                    disp(['Search area too small - increasing search area to ', num2str(sbnd_E_srch_act/meV), ', E_step=', num2str(E_step/meV), ', ', num2str(num_steps), ' steps']);

                else

                    k_t_mat_size    = size(k_t_mat, 2);
                    k_t_step_local  = k_t - k_t_mat(jj);
                    k_t_inter       = k_t - k_t_step_local/2;
                    k_t_mat         = [k_t_mat(1:jj), k_t_inter, k_t_mat(jj+1:k_t_mat_size)];

                    P               = polyfit(k_t_mat((jj-3):jj)/k_t_step, subband_E((jj-3):jj)/meV, 3);
                    E_predict       = polyval(P, k_t_mat(jj+1)/k_t_step)*meV;
                    dE              = E_predict - subband_E(jj);
                    sbnd_E_srch_act = sbnd_E_srch_old;
                    E_step          = E_step_old;

                    disp(['Inserting additional k_t=', num2str(k_t)]);
                end

                if (sbnd_E_srch_act>0.25*meV)
                    k_t_mat_size    = size(k_t_mat, 2);
                    k_t_step_local  = k_t - k_t_mat(jj);
                    k_t_inter       = k_t - k_t_step_local/2;
                    k_t_mat         = [k_t_mat(1:jj), k_t_inter, k_t_mat(jj+1:k_t_mat_size)];
                    P               = polyfit(k_t_mat((jj-3):jj)/k_t_step, subband_E((jj-3):jj)/meV, 3);
                    E_predict       = polyval(P, k_t_mat(jj+1)/k_t_step)*meV;
                    dE              = E_predict - subband_E(jj);
                    sbnd_E_srch_act = sbnd_E_srch_old;
                    E_step          = E_step_old;
                    disp(['Search area too large, inserting additional k_t=', num2str(k_t)]);
                end

                error = 0;
            end
        else                        % -- minimum found

            [f_h,f_l] = PlotWavefunction(E_result, k_t, 1, a);

            if ( ((f_l(N)^2 + f_h(N)^2)/(max(abs(f_l))^2+max(abs(f_h))^2)) < 1e-4*accuracy_modifier )
                subband_E(jj+1) = E_result;
                coeff(jj+1)          = a;

                if (jj<4)
                    dE          = (E_result - subband_E(jj))*(k_t_mat(jj+2)-k_t)/(k_t-k_t_mat(jj));
                elseif (jj<(size(k_t_mat,2)-1))
                    P           = polyfit(k_t_mat((jj-2):(jj+1))/k_t_step, subband_E((jj-2):(jj+1))/meV, 3);
                    E_predict   = polyval(P, k_t_mat(jj+2)/k_t_step)*meV;
                    dE          = E_predict - E_result;
                end

                for (cb_num=1:num_cb_subbands)
                    overlap1 = trapz(z, (wavcon(:,cb_num).*f_h'));
                    overlap2 = trapz(z, (wavcon(:,cb_num).*f_l'));

                    rtsTE(nr, cb_num,jj+1) = 0.5*(abs(overlap1^2) + (1/3)*abs(overlap2^2));
                    rtsTM(nr, cb_num,jj+1) =(2/3)*abs(overlap2^2);
                end  % -- end of for loop

                old_error = E_result - min_E - sbnd_E_srch_act;
                jj        = jj + 1;

                disp(['k_t=', num2str(k_t), ', E_acc=', num2str(E_acc_act/meV), ...
                    ', srch_dE=', num2str(sbnd_E_srch_act/meV), ' E_step=', num2str(E_step/meV) ...
                    ' E found ', num2str(E_result/meV), ' meV']);

                sbnd_E_srch_old  = sbnd_E_srch_act;
                sbnd_E_srch_act  = min([max(sbnd_E_srch_act/5, abs(5*old_error)), subband_E_search]);
                sbnd_E_srch_fnd  = sbnd_E_srch_act;

                if (sbnd_E_srch_act < sbnd_E_lim)
                    sbnd_E_srch_act = sbnd_E_lim;
                end

                sbnd_E_lim = max([sbnd_E_lim/10, abs(5*old_error), min(E_acc,E_acc_act*1e4)]);

                E_step_old = E_step;
                E_step     = sbnd_E_srch_act/10;
                E_step_fnd = E_step;

                E_acc_prev = E_acc_act;
                E_acc_act  = min(E_acc*1e5, E_acc_act*1.25);
                E_acc_fnd  = E_acc_act;

                accuracy_modifier = max(accuracy_modifier*0.75, 1);
                changed_num       = 0;

            else

                % Adjusting search parameters
                if (jj<3)
                    E_acc_act = E_acc_act/10;

                    if (E_acc_act<E_acc/20)
                        if (min_num==1)
                            min_num = 2;
                        else
                            min_num = 1;
                        end
                    end
                else
                    E_acc_act = E_acc_act*0.5;
                    disp(['E_acc changed to ', num2str(E_acc_act/meV)]);
                    accuracy_modifier = min(accuracy_modifier*2, 100);
                    disp(['accuracy_modifier=', num2str(accuracy_modifier)]);

                    if (E_acc_act<E_acc_prev/3)

                        k_t_mat_size   = size(k_t_mat, 2);
                        k_t_step_local = k_t - k_t_mat(jj);
                        k_t_inter      = k_t - k_t_step_local/2;
                        k_t_mat        = [k_t_mat(1:jj), k_t_inter, k_t_mat(jj+1:k_t_mat_size)];

                        if (jj>=3)
                            if (jj==3)
                                index = 1;
                            else
                                index = jj-3;
                            end
                            P         = polyfit(k_t_mat(index:jj)/k_t_step, subband_E(index:jj)/meV, 3);
                            E_predict = polyval(P, k_t_mat(jj+1)/k_t_step)*meV;
                        end

                        dE              = E_predict - subband_E(jj);
                        sbnd_E_srch_act = sbnd_E_srch_old;
                        E_step          = E_step_old;

                        disp(['Inserting additional k_t=', num2str(k_t_inter)]);
                    end

                    if (E_acc_act<E_result*eps*10)
                        if (~changed_num)
                            if (min_num==1)
                                min_num = 2;
                            else
                                min_num = 1;
                            end

                            E_step          = E_step_fnd;
                            sbnd_E_srch_act = sbnd_E_srch_fnd;
                            E_acc_act       = E_acc_fnd;
                            changed_num     = 1;
                        else
                            subband_E(jj+1) = E_result;
                            disp('Solution not found, assumed calculated value');

                            old_error = E_result - min_E - sbnd_E_srch_act;
                            jj        = jj + 1;

                            sbnd_E_srch_old = sbnd_E_srch_act;
                            sbnd_E_srch_act = min([max(sbnd_E_srch_act, abs(5*old_error)), subband_E_search]);

                            if (sbnd_E_srch_act<sbnd_E_lim)
                                sbnd_E_srch_act = sbnd_E_lim;
                            end

                            sbnd_E_lim = max([sbnd_E_lim/10, abs(5*old_error), min(E_acc,E_acc_act*1e4)]);

                            E_step_old = E_step;
                            E_step     = sbnd_E_srch_act/10;

                            E_acc_prev = E_acc_act;
                            E_acc_act  = E_acc_act*1e2;

                            [f_l, f_h] = PlotWavefunction(E_result, k_t, 1, a);

%                             figure;
%                             plot(f_h); hold; plot(f_l, 'r'); drawnow;

                            accuracy_modifier = min(accuracy_modifier*2, 100);

                        end
                    end
                end
            end
        end
    end  % -- end of while loop --

    t(nr) = toc;
    disp(['Elapsed time = ', num2str(t(nr)), ' sec']);

    % Plotting
    figure;
    plot(k_t_mat, subband_E/meV);
    title('E-k');
    xlabel('k_t (m^-^1)'); ylabel('E (meV)');
    drawnow;

    figure;
    plot(k_t_mat, squeeze(rtsTE(nr,1,1:size(k_t_mat,2))),cl(1));
    hold;
    for (cb_num=2:num_cb_subbands)
        plot(k_t_mat, squeeze(rtsTE(nr,cb_num,1:size(k_t_mat,2))),cl(cb_num));
    end    % -- end of for loop
    title('Transmition Strength - k_t');
    xlabel('Transmition Strength'); ylabel('k_t (m^-^1');

    figure;
    plot(k_t_mat, squeeze(rtsTM(nr,1,1:size(k_t_mat,2))),cl(1));
    for (cb_num=2:num_cb_subbands)
        plot(k_t_mat, squeeze(rtsTM(nr,cb_num,1:size(k_t_mat,2))),cl(cb_num));
    end    % -- end of for loop
    title('Transmition Strength - k_t');
    xlabel('Transmition Strength'); ylabel('k_t (m^-^1');
    
    eval(['save coeff' num2str(nr) ' coeff']);

end  % -- end of for loop --

%eval(['save rtsTE rtsTE']);
%eval(['save rtsTM rtsTM']);
